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Abstract 

In  this  paper,  an  analogy  between  the  mathematical 
modeling  of  transonic  potential  flow  and  the  flow  in  a 
cavitating  bearing  is  described.  Based  on  the  similarities, 
characteristics  of  the  cavitated  region  and  jump  conditions 
across  the  film  reformation  and  rupture  fronts  are  devel¬ 
oped  using  the  method  of  weak  solutions.  The  mathemati¬ 
cal  analogy  is  extended  by  utilizing  a  few  computational 
concepts  of  transonic  flow  to  numerically  model  the 
cavitating  bearing.  Methods  of  shock  fitting  and  shock 
capturing  are  discussed.  Various  procedures  used  in  tran¬ 
sonic  flow  computations  are  adapted  to  bearing  cavitation 
applications,  for  example,  type  differencing,  grid  trans¬ 
formation,  an  approximate  factorization  technique,  and 
Newton’s  iteration  method.  These  concepts  have  proved 
to  be  successful  and  have  vastly  improved  the  efficiency 
of  numerical  modeling  of  cavitated  bearings. 

Introduction 

Cavitation  in  fluid  film  bearings,  though  recognized  as 
early  as  1886  when  Reynolds  introduced  the  theory  of 


lubrication,  is  still  a  subject  which  draws  intense  debate 
as  to  the  nature  and  mechanism  of  the  phenomena.  Vari¬ 
ous  theories  and  conditions  for  cavitation  have  been  put 
forward.  However,  only  the  collective  works  of  Jakobsson 
and  Floberg  (1957)  and  Olsson  (1965),  now  known  as 
JFO  theory,  have  provided  insight  into  the  subject,  which 
is  both  consistent  with  mass  conservation  and  the  physics 
of  the  problem.  When  the  boundary  conditions  developed 
in  JFO  theory  are  applied  to  the  Reynolds  equation,  the 
extent  of  cavitated  regions  and  the  performance  of  bear¬ 
ings  can  be  predicted  more  precisely  than  by  any  existing 
method.  This  theory  has  yielded  results  which  are  in  good 
agreement  with  experimental  data. 

The  subject  of  gas  dynamics  gained  immense  research 
interest  around  the  turn  of  this  century.  This  effort  helped 
promote  the  development  of  supersonic  aircrafts  and  large 
thrust  rocket  nozzles.  Application  of  gas  dynamics  prin¬ 
ciples  include  turbine  flows,  gas  lasers,  aerodynamic 
windows,  missile  aerodynamics,  jet  engines  and  the  flow 
around  a  body  entering  the  atmosphere  (Emanuel,  1986). 
Mach  number  M,  a  nondimensional  parameter,  which  is 
the  ratio  of  the  flow  speed  to  the  local  speed  of  sound,  is 
the  indicative  index  as  to  whether  the  flow  is  subsonic 
(M<  1),  sonic  (M=  1 ),  or  supersonic  (M>  1).  It  is  also  a 
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measure  of  the  compressibility  of  the  flow.  In  a  converging- 
diverging  duct,  the  flow  can  range  from  subsonic  to 
supersonic  or  the  reverse.  The  converging  section  is  called 
a  nozzle  when  the  flow  is  subsonic  and  is  accelerating;  it 
is  called  a  diffuser  when  the  flow  is  supersonic  and  is 
decelerating.  The  diverging  section  is  called  a  diffuser 
when  the  flow  is  subsonic  and  is  decelerating  and  is  called 
a  nozzle  when  the  flow  is  supersonic  and  is 
accelerating.  When  subsonic  and  supersonic  flow  regimes 
exist,  the  flow  is  called  transonic.  It  is  also  possible  for  an 
internal  flow  to  be  totally  subsonic  or  supersonic  through 
the  nozzle.  However,  for  a  transonic  flow  to  exist,  a  duct 
with  a  throat  is  essential. 

Prior  to  1965,  computational  methods  were  rarely  used 
in  aerodynamic  analysis  and  importance  was  placed  on 
expensive  and  time  consuming  wind  tunnel  experiments. 
With  the  emergence  of  powerful  computers,  computa¬ 
tional  aerodynamics  lias  gieatly  auvanced  to  the  extent 
that  the  flow  pattern  past  entire  aircraft  in  different  flight 
regimes  can  be  predicted  (Jameson.  1987).  Such  rapid 
growth  in  computational  tcchniquescan  be  attributed  to  its 
direct  application  in  the  design  of  aircraft  and  space 
vehicles.  In  addition,  there  was  little  recourse  aside  from 
expensive  experimentation,  due  to  the  nonlinear  nature  of 
the  governing  equations  which  made  them  intractable  to 
analytical  modeling.  Developments  in  computational 
methods  applied  in  the  lubrication  area  have  been 
comparatively  slow  to  evolve.  In  fact,  to  date,  the  numeri¬ 
cal  algorithms  developed  by  Elrod  ( i 98 1 )  and  Kumar  and 
Booker  (1989)  are  the  only  effective  numerical  tools 
available  in  the  analysis  of  cavitation  in  bearings. 

Transonic  flow  theory  and  the  theory  of  lubrication  arc 
two  distinctly  different  fields  as  far  as  the  physical  phe¬ 
nomena  arc  concerned.  While  transonic  flow  theory  deals 
with  compressible  fluid  flow  near  sonic  velocity,  classical 
lubrication  theory  is  generally  concerned  with  the  flow  of 
a  highly  viscous  incompressible  fluid  with  a  Reynolds 
number  that  is  very  small.  However,  due  to  the  existence 
of  subsonic  and  supersonic  flow  regimes  in  a  converging- 
diverging  nozzle  and  the  cxistenccof  full  film  and  cavitatcd 
regions  in  a  bearing  with  a  converging-  diverging  wedge, 
there  exists  a  striking  resemblance  in  the  mathematical 
modeling  of  these  two  problems.  Such  an  analogy  can 
substantially  benefit  cither  field  by  suitably  incorporating 
the  advancements  from  one  field  to  the  other. 

In  this  paper,  a  mathematical  formulation  for  a  cavitatcd 
bearing  is  derived  and  compared  with  that  of  transonic 
potential  flow.  An  analogy  between  these  formulations 
are  developed.  The  analogy  is  utilized  to  employ  the 
method  of  characteristics  and  the  method  of  weak  solu¬ 
tions  from  transonic  flow  theory  to  determine  the 
characteristics  of  a  cavitatcd  region  and  the  jump  condi¬ 
tions  that  apply  across  both  film  reformation  and  rupture 
fronts  Tie  method  oi  determining  the  location  of  film 


reformation  using  shock  filling  and  shock  reformation 
techniques  are  discussed.  A  brief  discussion  of  several 
techniques  that  arc  widely  used  in  current  transonic  flow 
computation  is  provided.  These  techniques  have  already 
been  suitably  modified  and  incorporated  into  the  analysis 
of  cavitation  in  bearings. 

Mathematical  Modeling 


Cavitated  Bearing 


The  conservation  of  mass  flow,  within  the  clearance 
between  the  stationary  and  moving  surfaces  of  a  bearing 
can  be  written,  by  lumping  across  the  film  thickness,  as 


9ph 

"at" 


+  V  •  m  =  0 


0) 


In  the  converging  wedge  of  the  bearing,  the  film  thickness 
diminishes  and  the  pressure  is  developed.  In  this  region, 
the  mass  flux  mf  can  be  represented  by 


(2) 


The  first  term  on  the  right  side  is  the  mass  flow  due  to 
shear  (Couette  flow)  and  the  second  term  is  the  flow  due 
to  pressure  gradient  (Poiseuille  flow).  Somewhere  in  the 
diverging  wedge  of  the  bearing,  the  film  ruptures  and  a 
cavitation  region  is  formed  which  continues  until  the  film 
is  reformed  again.  In  this  cavitated  region,  the  pressure 
remains  essentially  constant  at  the  cavitation  pressure  and 
the  mass  flows  across  this  region  due  only  to  shear  along 
the  film  striatiens.  In  the  cavitated  region,  the  film  occu¬ 
pies  only  a  portion  of  the  volume,  the  remaining  portion 
being  filled  by  air,  gas,  or  vapor.  Mass  flow  in  this  region 
is  given  by 


mc  = 


(3) 


where  0  is  the  partial  film  content  in  the  cavitatcd  region. 
Although  the  film  consists  of  incompressible  fluid,  if  it  is 
assumed  that  the  density  of  the  film  varies  due  to  the 
applied  pressure,  then  the  variable  0  can  be  provided  with 
a  dual  interpretation 


0  = 


p/pc,  in  the  full  film  region  where  0  >  1 
Vf  /  V, ,  in  the  cavitated  region  where  0  <  1 


where  Vf  and  V,  arc  the  volume  occupied  by  the  film  and 
the  total  volume,  respectively.  In  the  full  film  region,  the 
variation  of  density  will  be  governed  by  the  bulk  modulus 
of  the  liquid,  that  is. 


P  = 


pf*R  =  0^E 
^9p  90 


(4) 
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Equation  (4)  also  enables  one  to  represent  the  pressure  in 
terms  of  the  density  (or  0),  and  in  essence,  acts  as  the 
equation  of  state  of  the  lubricant 
Mass  flow  through  theentire  bearing  can  then  be  written 
as 

m  =  -pch0  -  g V0  (5) 

2  K  6  12)i 

where  g  is  a  switch  function,  which  is  introduced  to 
remove  the  flow  due  to  pressure  gradient  within  the 
cavitated  region  and  is  defined  by 


8  = 


1  when  9  >  1 
0  when  a  <  1 


Fora  finite  bearing,  the  flow  due  to  shear  occurs  only  in  the 
circumferential  direction  while  the  flow  due  to  pressure 
gradient  is  present  in  both  the  circumferential  and  axial 
directions.  Hence,  the  two  dimensional  form  of  equation  ( 1 ) 
can  be  written  as 

3pc0h  3  f  pchU9  pc(3h3  30 ' 

3t  3x  2  12p  *  3x^ 
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PcPh3  ae' 

12p  3z^ 
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Equation  (6),  which  can  be  used  to  describe  the  mass  flow 
through  the  entire  bearing,  was  developed  by  Elrod  and 
Adams  (1974)  and  termed  a  ‘universal  equation’.  In  the 
full  film  region,  equation  (6)  may  be  written  as 


3E  U  3E 
3t  2  3x 


where  E  =  0h  and  K=  -|ih3g/12|i.  Equation  (7)  is  an 
elliptic  partial  differential  equation.  The  form  of  equa¬ 
tion  (6)  in  the  cavitated  region  is  given  by 


3E  V  3E 
3t  +  2  3x 


=  0 


(8) 


and  is  a  hyperbolic  partial  differential  equation.  This  may 
be  easily  demonstrated  by  differentiating  with  respect  to 
t  to  get 

32E  U  32E  32E  U2  32E  „ 

3t2  2  3x3t  3t2  4  3x2 

Equation  (9)  is  a  canonical  form  of  the  wave  equation.  The 
characteristic  form  of  this  equation  has  two  real  roots, 
±U/2. 

In  the  full  film  region,  pressure  increases  toa  maximum 
value  and  then  gradually  decreases  until  the  pressure  and 
its  derivative  simultaneously  vanish,  at  a  location  where 
the  film  ruptures.  The  air/gas/vapor  strips  in  the  cavitated 
region  begin  with  a  pointed  shape.  When  the  film  reforms, 
depending  upon  the  upstream  conditions,  the  reformation 
occurs  abruptly.  This  is  because  the  cavitated  region  is  not 
able  to  signal  the  impending  conditions  to  the  upstream 
flow.  Figure  1  represents  typical  profiles  of  pressure  and 
fractional  film  content  for  a  submerged  journal  bearing,  at 
a  particular  axial  location.  The  abrupt  changes  in  0  and 
the  pressure  gradient  at  the  reformation  boundary  and  the 
gradual  changes  at  the  rupture  boundary  can  be  clearly 


seen.  When  the  effects  of  cavitation  arc  not  formally 
considered,  it  is  usually  assumed  that  the  film  ruptures 
at  the  minimum  film  thickness  and  reforms  at  the  maxi¬ 
mum  film  thickness.  When  the  cavitation  boundaries  are 
determined,  it  is  found  that  the  film  extends  slightly 
beyond  the  minimum  film  thickness  into  the  diverging 
portion  of  the  bearing  and,  depending  upon  the  lubricant 
supply  conditions,  the  film  reformation  can  occur  at  or 
around  the  feed  groove. 

Transonic  Flow 


The  flow  of  a  compressible  fluid  in  thermodynamic 
equilibrium  is  governed  by  the  Navicr-Stokes  equations. 
Fora  two-dimensional  flow  these  equations  can  be  written 
in  vector  form  as 


3w  +  3f  3g  _  3R  +  3S 
3t  3x  3y  3x  3y 


(10) 


where  w  is  the  vector  of  dependent  variables:  density, 
Cartesian  velocity  components,  and  total  energy;  f  and  g 
are  the  convective  flux  vectors;  and  R  and  S  are  the 
viscous  flux  vectors.  Because  the  full  Navicr-Stokcs 
equations  arequite  complex, approximationsare  generally 
made.  One  such  simplification  consists  of  assuming  no 
viscous  dissipation  and  that  flow  is  irrotational.  Conse¬ 
quently,  equation  (10)  can  be  written  in  a  quasi-linear 
form  in  terms  of  the  velocity  potential  0. 
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It  should  be  noticed  that,  in  a  subsonic  flow  regime 
(u,v  <  c),  the  coefficients  of  the  second  order  terms  will  be 
positive  and,  in  a  supersonic  flow  regime  (u,v  >  c),  the 
coefficients  become  negative.  This  variation  results  in  the 
equation  being  of  the  elliptic  type  with  two  imaginary 
roots  in  subsonic  regions  and  of  the  hyperbolic  type  with 
two  real  roots  in  supersonic  regions. 

When  a  subsonic  flow  slows  down,  it  does  so  gradually. 
On  the  other  hand,  a  supersonic  How,  w'hich  can  also 
decelerate  gradually,  normally  slows  down  abruptly. 
Because  the  fluid  in  a  supersonic  flow  is  unable  to  signal 
the  upstream  flow  of  any  flow  oi  geometric  changes.  This 
is  a  typical  characteristic  of  phenomena  governed  by 
hyperbolic-type  equations.  The  abrupt  changes  lead  to 
discontinuities  in  the  flow  which  are  known  as  shock 
waves  and  arc  a  distinct  feature  of  a  supersonic  flow  in 
establishing  the  overall  nature  of  the  flow  field.  Of  course, 
an  internal  flow  can  also  emerge  as  supersonic  without  the 
presence  of  any  shock  if  the  outlet  conditions  permit. 
Since  viscous  effects  arc  neglected  in  the  potential  equa¬ 
tion,  for  internal  flow,  sonic  conditions  occur  at  the  throat 
of  the  ducL  For  real  fluids,  the  sonic  line  extends  slightly 
downstream  of  the  throat  into  the  diverging  portion  of  the 
duct. 

Computational  techniques  for  potential  flow  have  been 
extensively  developed,  since,  extremely  inexpensive 
computation  is  achieved  with  this  formulation.  Moreover, 
shock  capturing  and  convergence  acceleration  techniques 
developed  for  potential  flow  have  been  found  to  be  trans¬ 
ferable  to  more  complex  models  using  Euler  equations. 


Film  reformation 
(with  cavitation) 


Film  reformation 


Figure  2. — Transonic  flow  and  cavitatcd  bearing. 
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Analogy 

The  similarities  between  internal  transonic  flow  and 
cavitated  bearing  modeling  are  evident  from  the  previous 
two  sections.  Figure  2  illustrates  these  similarities.  Both 
phenomena  are  governed  by  similar  mathematical  formu¬ 
lations  in  the  different  regions  and  both  have  an  embedded 
hyperbolic  region  within  the  elliptic  region.  The  subsonic 
portion  of  the  flow  (M  <  1)  is  analogous  to  the  full  film 
region  (0  >  1).  The  sonic  line  (M  =  1)  is  analogous  to  film 
rupture  (0  =  1).  The  supersonic  portion  of  the  flow  (M  >  1) 
is  analogous  to  the  full  film  region  (0  <  1).  Also,  a  shock 
and  film  reformation  have  similar  characteristics.  It  should 
be  recognized  that  the  Mach  number  and  the  inverse  of 
fractional  film  content  have  similarities.  The  sonic  line/ 
the  film  rupture  locations  are  strongly  influenced  by  the 
geometry  of  the  flow,  while  the  shock  wave/film  reforma¬ 
tion  locations  are  due  to  the  upstream  conditions.  The  flow 
can  also  be  fully  elliptic  or  hyperbolic  in  both  cases, 
although  a  completely  cavitated  bearing  has  no  physical 
significance.  Similar  to  compressible  flows  involving 
shocks,  determination  of  the  film  reformation  boundary  is 
a  difficult  task. 

Although,  it  is  seen  that  transonic  flow  and  flow  in  a 
cavitated  bearing  have  similarities,  it  should  also  be  noted 
that  they  also  differ  in  several  respects.  The  essential 
differences  between  these  two  models  are  the  following: 

( 1 )  In  transonic  flow,  the  type  of  the  equation  is  changed 
due  to  the  change  in  the  sign  of  the  coefficient  of  the 
second  order  terms;  in  the  case  of  a  cavitated  bearing,  the 
second  order  terms  are  totally  lost  in  the  hyperbolic  region 
resulting  in  the  reduction  of  the  order  of  the  governing 
equation.  This  sometimes  causes  oscillations  at  the 
boundary  locations. 

(2)  The  entire  flow  is  compressible  in  a  transonic  flow; 
but,  in  a  cavitated  bearing,  the  full  film  region  flow 
although  really  incompressible  is  taken  to  be  compress¬ 
ible  and  the  cavitated  region  flow  is  of  the  compressible 
type. 

(3)  The  unknowns  in  an  irrotational  transonic  flow  are 
density  and  the  potential  function  which  are  dependent  on 
each  other.  On  the  other  hand,  for  a  cavitated  bearing  the 
only  unknown  is  density  (or  0). 

(4)  Row  in  both  Cartesian  directions  can  occur  in 
transonic  flow,  although  the  resultant  velocity  can  be 
made  to  align  with  one  of  the  coordinate  axes  by  Jameson ’s 
rotated  difference  scheme  (1974).  In  bearings,  generally 
the  only  motion  is  a  direct  result  of  the  journal  rotation  and 
the  velocity  vector  is  normally  taken  to  coincide  with  a 
coordinate  axis. 


Having  pointed  out  the  basic  analogy  between  the  two 
models,  extension  of  method  of  characteristics  to  deter¬ 
mine  the  path  of  the  disturbances  in  the  cavitated  region 
and  the  method  of  weak  solutions  to  determine  the  jump 
conditions  across  the  discontinuities  can  be  developed. 
Also,  determination  of  the  location  of  film  reformation 
using  shock  fitting  and  shock  capturing  methods  will  be 
discussed.  Since,  our  intention  is  to  develop  correspond¬ 
ing  expressions  for  a  cavitated  bearing  based  on  the 
analogy,  the  development  details  for  transonic  flow  are 
not  presented  in  detail,  here.  Interested  readers  are  refered 
to,  for  example,  Anderson  et  al.  (1984). 

Method  of  Characteristics 

Hyperbolic  equations  have  certain  lines  (or  surfaces) 
which  indicate  the  zones  of  influence  and  zones  of 
dependence.  The  information  about  the  flow  is  signalled 
along  these  lines  which  are  called  characteristics.  This 
property  is  used  to  determine  the  value  of  the  variable  at 
a  particular  location  in  a  hyperbolic  region  from  the 
known  value  of  the  variable  at  a  downstream  location. 
This  method  of  solving  hyperbolic  equations  is  the  Method 
of  Characteristics  (MOC). 

Supersonic  flow.  -  Assuming  the  free  stream  is  aligned 
with  the  x  axis,  equation  (1 1)  can  be  written  as 

(l  •  ML]4>XX  +  <|>yy  =  0  (12) 

where  M„  is  the  freestream  Mach  number.  In  order  to 
determine  the  characteristic  direction,  equation  (12)  is 
written  in  terms  of  the  Cartesian  velocity  components 
along  an  arbitrary  smooth  curve  C,  and  the  determinant  of 
the  coefficients  is  set  to  zero  along  the  curve.  This  will 
result  in  the  differential  equations  for  the  characteristics. 
For  this  case, 


dy  _  ±1 

dx  B 

B2  =  (mI  - 1) 


(13) 


For  a  constant  B,  the  equations  describing  the  character¬ 
istics  are  obtained  by  integrating  equation  (13).  This 
results  in  the  following  form: 


4  =  x  -  By 
T)  =  x  +  By 


(14) 
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where  %  and  q  are  coordinates  along  the  characteristic 
line. 

Cavitated  region.  -  Within  the  cavitated  region,  the 
governing  hyperbolic  equation  is  written  as 


This  eliminates  the  need  for  the  solution  to  be  differen¬ 
tiable  across  the  discontinuity.  The  mathematical  theory 
of  weak  solutions  for  hyperbolic  equations  is  a  relatively 
recent  development  and  may  be  utilized  to  determine  the 
jump  conditions  across  a  discontinuity  in  a  flow. 


3E  _U  3F 

3t  2  3x 


=  0 


(15) 


Along  a  curve  x  =  x(t)  in  the  x  - 1  plane,  E  =  E(x).  For  a 
particular  curve  x  =  xc(t),  let  dE  =  0. 

Thus, 


Transonic  flow.  -  Although  the  steady  state  formula¬ 
tions  exhibit  elliptic  and  hyperbolic  type  equations  at 
different  parts  of  the  flow  field,  addition  of  an  unsteady 
term  results  in  hyperbolic  type  of  equation.  Consider  a 
one-dimensional,  scalar,  nonlinear,  hyperbolic  partial 
differential  equation 


dE  = 


3E 
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Using  equation  (15),  results  in 
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The  solution  to  the  first  of  equation  (18)  determines  the 
equation  of  the  characteristics  and  the  second  one  reveals 
the  parameter  that  is  constant  along  the  characteristics. 
They  are 


k  =  xc  - 


E  =  Oh  =  constant 


(19) 


Since  h  is  known,  the  value  of  0  at  any  point  in  the 
cavitated  region  can  be  obtained  from  a  point  with  a 
known  value  of  0  by  tracing  backwards  along  the  char¬ 
acteristic  line.  Olsson  (1965)  discussed  this  characteristics 
approach  in  his  treatise  on  dynamically  loaded  bearings. 


Method  of  Weak  Solutions 


A  genuine  solution  of  a  hyperbolic  differential  equation 
is  one  in  which  the  dependent  variable  is  continuous  but 
discontinuities  in  its  derivatives  may  occur.  Alternately ,  a 
weak  solution  is  genuine  except  along  a  surface  across 
which  the  dependent  variable  is  discontinuous.  Only  the 
scalar  or  vector  form  of  first  order  and  hyperbolic  second 
order  partial  differential  equations  possess  weak  solu¬ 
tions.  Since  the  dependent  variable  is  not  continuous 
across  the  discontinuity,  an  integral  formulation  is  used. 


where  u  and  F(u)  are  unknown  variables.  This  can  be 
rewritten  as 


3u 

aT 


A 

+  A37 


■  0 


(21) 


where  A  =  A(u)  =  dF/du  is  called  the  Jacobian.  Now,  if 
w(x  ,t)  is  an  arbitrary  test  function  which  is  continuous  and 
has  a  continuous  first  derivative  but  vanishes  on  the 
boundary  and  outside  of  an  arbitrary  domain  D  inthe(x,t) 
plane,  then 


JJD(f 

When  both  u  and  F  are  continuous  and  have  continuous 
first  derivatives,  it  can  be  shown  that 


JJD(uwt  +  Fwx)dxdt  =  0  (23) 

Functions  u(x,t)  which  satisfy  equation  (23)  for  all  lest 
functions  w  are  called  weak  solutions  of  equation  (21).  If 
the  domain  D  contains  a  moving  curve  x(x,t),  across  which 
u  is  discontinuous  as  shown  in  figurc3,  equation  (23)  may 
be  integrated  by  parts,  utilizing  equation  (22),  to  get 

j^([u]cosai  +  [F]cosu2)ds  =  0  (24) 

where  [  ]  denotes  a  jump  across  the  discontinuity  and 
cos  a,  and  cos  are  the  direction  cosines  normal  to  the 
discontinuity  x  (x,t).  Since  the  integrand  must  vanish  for 
all  w,  the  condition  for  u  to  be  a  weak  solution  of 
equation  (21)  will  therefore  be 


[ujcosai  +  [F]cosot2  =  0  (25) 

Equation  (25)  determines  the  jump  in  the  value  of  u 
across  the  discontinuity.  This  analysis  is  also  valid  for  a 
system  of  equations  in  which  case,  u  and  F  are  vectors. 
The  jump  conditions, 
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Figure  3. — An  arbitrary  domain  with  a  discontinuity. 
(26) 


where  dx/dl  is  the  velocity  of  the  discontiuity,  are 
popularly  known  as  the  Rankine-Hugoniot  equations.  For 
example,  if  u  and  F  are  vectors  defined  as 


where  A  =  A(E)=dm/dE.  To  proceed  in  the  same  manner 
as  that  used  for  the  transonic  flow  illustration,  it  can  be 
shown  that  for  E  to  be  a  weak  solution  of  equation  (28), 
the  jump  condition  to  be  satisfied  at  any  discontinuity  is 


u 


p,pu  pe  +  Py 


[E]  cos  oti  +  |m]  cos  a2  =  0 


(30) 
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then  the  Rankine-Hugoniot  equations  can  be  written  as 

ud[p]  =  [p“] 
ud[pu]  =  [p«2  +  p] 

Ud[pe  +  pu2  /  2]  =  [u(pe  +  (pu2/2)  +  p)j 

where  Ud  is  the  velocity  of  discontinuity  and  e  is  the 
internal  energy.  The  first  two  equations  are  called 
mechanical  jump  conditions. 

Cavitated  bearing.  -  Consider  the  one-dimensional 
version  of  equation  (6),  which  can  be  written  as 


At  any  time  t] ,  the  velocity  of  propagation  of  the  discon¬ 
tinuity  is  determined  by  the  tangent  of  the  curve  as  shown 
in  figure  4.  The  direction  cosines  are 


cosotj  = 


therefore, 


1  + 


-and  cosa2  =  -■ 
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dx 
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dxA=t, 
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cos  ct)  _  dx 
cos  a2  dt 

Hence,  equation  (30)  can  be  written  in  terms  of  the 
primitive  variables  as 

Ml  -  Mr]£  -  JfK  -  Mr] 

!*Ol  ■ 


12p 


=  0  (32) 


where  m  is  the  mass  flux  as  defined  by  equation  (5). 
Equation  (28)  can  be  written  in  a  form  similar  to  that  of 
equation  (21),  that  is, 


where  L  and  R  are  the  conditions  at  the  left  and  right  sides 
of  the  discontinuity,  respectively.  In  other  words, 


8E  *  aE  n 

¥  +  asT  =0 


(29) 


(»]fj  -  H 

is  an  equivalent  Rankine-Hugoniot  equation. 


(33) 
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Figure  4. — Path  of  a  discontinuity. 


First  consider  the  discontinuity  due  to  Film  reformation. 
If  it  is  assumed  that  the  film  thickness  variation  across  the 
discontinuity  is  negligible  and  that  0R=  1  since  full  film 
is  formed  at  the  right  side  of  the  discontinuity,  equa¬ 
tion  (32)  for  these  conditions  can  be  written  as 


(eL  -  0$ 


t(0l  -  - 


=  o 


(34) 


If  the  reformation  front  is  not  moving  with  respect  to  time 
(analogous  to  a  stationary  shock),  then  the  conditions  for 
film  reformation  will  be 


— (0L  -  1)  =  JL 

2V  L  ’  12 p 


.2  30 

h  £ 


(35) 


Since  (30/3x)  >  0,  obviously,  0L  <  1. 

On  the  other  hand,  if  the  reformation  front  is  in 
motion,  equation  (34)  can  be  written  as  follows 


Since  (30/3x)  >  0,  the  right  side  of  equation  (36)  is  less 
than  or  equal  to  zero.  Also,  since  0L<1,  there  can  be  two 
conditions 


(i) dx/dt  >  U/2 

For  this  condition,  the  left  side  of  equation  (36)  is  greater 
than  or  equal  to  zero.  Hence,  the  only  condition  that  will 
satisfy  the  equality  is  that  both  sides  must  be  zero,  that  is, 
(30/3x)r  =  o  and  0^  =  1.  This  is  the  classical  Reynolds 
boundary  condition  which  is  generally  applied  to  deter¬ 
mine  film  rupture. 


(ii)  dx/dt  <  II / 2 

Now,  the  left  side  of  equation  (36)  is  also  less  than  or 
equal  to  zero.  Hence,  qL  cannot  be  uniquely  determined 
from  equation  (36);  it  must  be  obtained  using  the  charac¬ 
teristics.  Equation  (36)  is  treated  as  a  differential  equation 
of  the  following  form  to  determine  the  velocity  of  the 
front: 

*  =  !L.± (37) 

dt  2  12p  (1  -  0L) V  3xJr  P  ; 

Although  in  transonic  flow,  the  sonic  line  is  not  treated 
as  a  discontinuity;  in  the  case  of  a  dynamically  loaded, 
cavitated  bearing,  there  could  be  a  discontinuity  of  0  at 
the  film  rupture,  when  dx/dt>U/2  and  there  will  bea  jump 
in  its  gradient.  The  jump  conditions  for  this  case  can  be 
developed  in  the  same  manner  as  previously  described. 
The  resulting  expressions  will  be 

(i)  dx/dt  S  U/2  (96/dx)L  -  Oand  e„  -  1 


(ii)  dx/dt  >  U/2 


dx 

dt 


H  +  JL  1  fh2*0 

2  12p  (e,  - 1)1  3xJt 


(38) 


These  jump  conditions  are  exactly  the  same  as  those 
determined  by  Olsson  (1965).  However,  in  that  reference 
the  conditions  were  derived  using  a  mass  balance  across  a 
fluid  volume  containing  the  discontinuity. 


Computational  Treatment  of  Discontinuity 


Transonic  flow.  -  The  numerical  computation  of 
supersonic  flow  is  complicated  due  to  the  presence  of 
shock  waves,  across  which  the  dependent  variables  and 
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their  derivatives  may  be  discontinuous.  Two  types  of 
numerical  techniques  have  been  developed  to  analyze 
such  flow  fields  and  are  known  as  shock  fitting  and  shock 
capturing  techniques. 

Shock  fitting  technique.  -  This  technique  attempts  to 
locate  any  discontinuities  and  treats  them  as  boundaries 
between  the  regions  of  the  flow  field  where  regular  solu¬ 
tions  are  applicable.  Shock  fitting  is  achieved  by  satisfying 
the  Rankine-Hugoniot  equations  across  the  discontinuities 
while  simultaneously  ensuring  that  the  solution  on  the 
downstream  side  of  any  shock  is  compatible  with  the  rest 
of  the  flow  field.  The  movement  of  the  shock  wave  is 
obtained  as  a  part  of  the  solution.  The  flow  field  down¬ 
stream  of  each  shock  can  be  determined  from  freestream 
conditions.  If  the  upstream  conditions,  initial  shock  slope 
and  velocity  arc  known,  the  shock  acceleration  and  post 
shock  pressure  can  be  determined  by  combining  Rankine- 
Hugoniot  equations  with  die  compatible  equation.  This 
technique  is  most  convenient  for  governing  equations 
written  in  nonconservative  form. 

Several  approaches  have  been  devised  to  fit  shocks 
(Moretti,  1974).  The  flow  field  is  either  partitioned  by 
aligning  any  shock  waves  with  grid  lines  or  the 
discontinuities  are  treated  explicitly,  but  not  as  bound¬ 
aries,  in  the  computation. 

Shock  capturing  technique.  -  Unlike  the  shock  fitting 
technique,  with  this  method,  the  discontinuities  are  pre¬ 
dicted  as  a  part  of  the  solution  without  the  requirement  of 
any  special  treatment.  By  casting  Euler  equations  in 
conservation-law  form,  the  weak  solutions  and  jump 
conditions  are  built  in.  The  conservative  form  of  govern¬ 
ing  equations  and  the  discretization  automatically  allow 
prediction  of  the  shock  wave  speed  and  the  strength 
(Lax,  1954). 

Because  of  the  simplicity  in  approach,  this  technique  is 
most  popular  in  the  computation  of  flows  with  shocks. 
This  technique  also  has  several  variants  which  include 
flux  splitting  and  split  coefficient  matrix  methods.  Shock 
waves  predicted  with  this  technique  can  be  smeared  over 
several  grid  spaces  and  the  application  of  surface  bound¬ 
ary  conditions  can  be  difficult.  Due  to  the  wavelike  nature 
of  hyperbolic  equations,  boundary  errors  are  propagated 
into  the  flow  field  which  results  in  instability  in  the 
computation.  In  general,  shock  capturing  techniques  arc 
applied  to  predict  internal  shocks  while  boundary  shocks 
are  fitted. 

Cavitated  bearing.  -  Similar  to  the  previously  described 
approaches  for  transonic  flow,  discontinuity  fitting  and 
capturing  techniques  can  be  applied  to  problems  with 
cavitating  regions.  Determining  the  film  reformation 


boundary  is  more  difficult  than  determining  the  film 
rupture  boundary  due  to  sudden  changes  in  the  flow 
variables  across  the  front. 

If  the  initial  location  and  slope  of  the  boundaries  are  not 
known,  they  can  be  determined  by  employing  a  trial  and 
error  method  (the  discontinuity  fitting  method).  The  loca¬ 
tion  of  the  boundaries  can  be  assumed  and  the  flow  field 
on  either  side  of  the  discontinuity  can  be  determined.  The 
equivalent  Rankine-Hugoniot  equation  can  then  be  ap¬ 
plied  across  the  boundary  to  verify  the  assumption.  If  the 
initial  location  and  slope  arc  known,  then  the  governing 
equations  coupled  with  the  equivalent  Rankine-Hugoniot 
equation  can  be  solved  to  determine  the  new  velocity  and 
the  new  locations  of  the  discontinuities. 

The  algorithm  introduced  by  Elrod  (1981)  is  essentially 
a  discontinuity  capturing  technique.  By  combining 
the  governing  equations  for  the  full  film  and  cavitated 
regions  and  conserving  mass  flow  through  the  entire 
bearing,  the  ‘  universal  equation’  is  cast  mtoa  conservation- 
law  form.  Hence,  the  discontinuities  can  be  predicted  as  a 
part  of  the  solution.  This  method  is  simple  to  implement 
and  docs  not  require  any  knowledge  about  the  location  of 
the  discontinuities.  The  boundaries  are  predicted  very 
effectively. 

Concepts  from  Transonic  Flow 
Computation 

The  authors  have  utilized  a  few  transonic  flow  compu¬ 
tational  concepts  in  the  analysis  of  cavitated  bearings.  The 
following  is  a  brief  discussion  of  this  work. 

Computational  Algorithm 

When  the  potential  flow  equation  was  used  in  transonic 
flow  computation,  difficulty  was  encountered  due  to  the 
reversal  of  the  velocity  vector  in  the  supersonic  flow 
regime.  Murman  and  Cole  (1971),  in  a  landmark  paper, 
demonstrated  a  simple  way  to  obtain  a  meaningful  solu¬ 
tion  by  proposing  the  use  of  central  differencing  in  the 
subsonic  region  (elliptic)  and  one  sided  upwind 
differencing  in  the  supersonic  region  (hyperbolic).  Jameson 
(1975)  created  a  type  differencing  scheme  by  introducing 
an  artificial  viscosity  term  into  the  governing  equation. 
This  enables  automatic  switching  of  the  form  of 
differencing  as  required  within  different  regions  of  the 
flow.  Also,  with  the  explicit  addition  of  an  artificial 
viscosity  term,  the  conservation  form  of  the  equation  was 
preserved. 

In  the  analysis  of  cavitation  in  bearings,  Elrod  (1981) 
modified  the  originally  proposed  algorithm  by  Elrod  and 
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Adams  (1974),  by  incorporating  the  idea  of  changing  the 
form  of  differencing  in  the  full  film  and  cavitated  regions. 
However,  this  effect  was  achieved  by  ‘trial  and  error’  and 
the  algorithm  was  empirically  developed.  A  type 
differencing  scheme  was  developed  by  Vijayaraghavan 
and  Keith  (1989),  by  introducing  an  artificial  viscosity 
term  (much  like  Jameson’s)  into  the  governing  equation 
which  in  turn  permitted  the  algorithm  to  be  mathemati¬ 
cally  derived.  In  addition,  with  this  modification,  the 
discretization  was  also  performed  in  conservative  form. 
The  predictions  using  this  modified  algorithm  were  found 
to  compare  well  with  the  predictions  using  the  Elrod 
algorithm  and  with  experimental  data.  Hence,  the  modi¬ 
fied  cavitation  algorithm,  with  this  firmer  base,  was  thought 
to  offer  greater  potential  for  further  improvement 

Grid  Control 

Numerical  grid  generation  is  a  fairly  comm  ,n  tool  to 
model  arbitrarily  shaped  regions  in  computational  fluid 
dynamics.  This  is  basically  a  procedure  to  distribute  in  an 
orderly  manner  the  grid  points  in  the  physical  field  in  such 
a  way  that  efficient  communication  between  the  points 
and  all  the  physical  phenomena  on  the  entire  continuous 
field  is  represented  with  sufficient  accuracy  (Thompson, 
1984).  Also,  the  region  in  the  immediate  vicinity  of  solid 
surfaces  are  dominant  in  determining  the  character  of  the 
flow  due  the  large  gradients  prevailing  in  this  region. 
Accurate  prediction  of  flow  variables  in  this  region  is 
important  since  the  final  values  of  the  variables  strongly 
depends  on  this  boundary  prediction.  When  such  high 
gradient  regions  are  not  known  d  priori,  a  dynamically 
adaptive  grid  system  can  be  an  effective  tool.  This  is  an 
active  area  in  grid  generation  research.  By  dynamically 
readjusting  the  grid  distribution  as  the  solution  proceeds, 
high  resolution  solutions  can  be  obtained  with  fewer  grid 
points. 

In  the  case  of  a  cavitated  bearing,  film  rupture  and 
reformation  locations  arc  not  known  &  priori.  Also,  accurate 
prediction  of  the  pressure  distribution  is  the  primary 
requirement  in  the  full  film  region,  since,  all  the  perform¬ 
ance  parameters  depend  upon  the  pressure  profile.  Hence, 
closely  placed  grid  points  around  the  cavitation  bound¬ 
aries  and  more  grid  points  in  the  high  pressure  gradient 
regions  should  lead  to  a  more  accurate  prediction  of  the 
pressure  profile.  In  the  case  of  a  bearing  with  a  misaligned 
journal,  the  maximum  pressure  location  in  the  axial  direc¬ 
tion  is  shifted  towards  the  edge  of  the  bearing  and  to 
correctly  predict  the  pressure  distribution  in  this  region, 
closely  spaced  grid  points  are  required.  With  a  uniformly 
fine  grid  arrangement,  this  will  result  in  a  large  number  of 
grid  points  located  within  regions  where  such  a  small  grid 


spacing  does  not  necessarily  provide  any  significant 
improvement  in  the  accuracy.  With  a  grid  adaption  tech¬ 
nique  all  these  effects  can  be  achieved  by  moving  the  grid 
points  and  selectively  loca_  ig  them  in  such  a  way  that 
accurate  solutions  can  be  obtained  with  fewer  grid  points. 
The  concept  of  gnd  generation/transformation  and  grid 
adaption  techniques  have  been  incorporated  into  the 
modified  cavitation  algorithm  by  Vijayaraghavan  and 
Keith  (1990(a)). 

The  grid  adaption  procedure  can  be  tailored  to  the 
problem  being  solved.  However,  the  procedure  envisaged 
was  to  perform  initial  computations  which  locate  the  film 
rupture  and  reformation  fronts,  to  distribute  closely  spaced 
grid  arcind  them,  and  then  to  rearrange  the  grid  in  the  full 
film  region  according  to  the  pressure  gradient.  In  the  case 
of  a  misaligned  journal,  when  the  degree  of  i  •  salignment 
is  large,  grid  adaption  in  the  axial  direction  is  applied  to 
cluster  the  grid  around  the  maximum  pressure  location 
(Vijayaraghavan  and  Keith,  1 990(b)>.  By  aligning  the  grid 
with  discontinuity  s,  the  flow  field  is  divided  into  zones  of 
full  film  and  cavitated  regions.  This  enables  the  jump 
conditions  to  be  applied  effectively.  In  addition,  in  the 
case  of  a  two-dimensional  time  asymptotic  solution,  by 
aligning  the  discontinuity  along  one  coordinate  direction, 
according  to  the  equivalent  Rankinc-Hugoniot  equation, 
the  mass  flow  along  the  other  coordinate  direction  is 
continuous  across  the  discontinuity.  This  method  of  grid 
adaption  combines  features  of  both  shock  fitting  and 
shock  capturing  techniques. 

The  predicted  performance  of  the  bearings  using  the 
adaptive  grid  and  conventional  orthogonal  grid  arrange¬ 
ments  were  found  to  be  comparable.  The  results  obtained 
thus  far  demonstrate  the  usefulness  of  these  techniques  in 
the  analysis  of  V.tring  problems.  The  transformation  of 
the  governing  equation  and  numerical  differencing  in  a 
nonorthogonal  coordinate  svsiem  could  be  confidently 
applied,  primarily  due  to  the  mathematical  base  provided 
in  the  modified  algorithm. 

Solution  Procedure 

In  the  case  of  transonic  flow,  for  steady  problems,  the 
converged  solution  obtained  by  using  relaxation  methods 
generally  requires  a  very  large  number  of  iterations  due  to 
the  slow  convergence  rate.  The  two  most  effective  solu¬ 
tion  acceleration  techniques  for  rapid  convergence  in  the 
transonic  flow  computations  are  approximate  factoriza¬ 
tion  of  the  difference  operators  and  the  use  of  multiple 
grids  (Jameson,  1987).  For  genuine  unsteady  transonic 
flow  problems,  Newton  iteration  tcchniqucscanbeapplicd 
to  the  governing  equation,  at  every  time  step,  to  obtain 
time  accurate  solutions  (Shankar  ct  al,  1985). 
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In  the  analysis  of  cavitatcd  bearings,  relaxation  meth¬ 
ods  are  found  to  also  require  a  large  number  of  iterations 
to  obtain  a  converged  steady  state  solution.  Woods  and 
Brewe  (1989),  in  the  analysis  of  a  dynamically  loaded, 
submerged  journal  bearing,  incorporated  a  multigrid  tech¬ 
nique  into  the  Elrod  algorithm  and  obtained  considerable 
savings  of  computer  lime,  h  an  unsteady  analysis  of 
cavitatcd  bearings,  particularly  uynamically  loaded  bear¬ 
ings,  time  accurate  solutions  are  very  important.  In  such 
cases,  the  accuracy  of  the  solution  can  be  improved  by 
adding  a  Newton  iteration  technique.  The  Newton  itera¬ 
tion  scheme  and  approximate  factorization  techniques 
were  developed  by  Vijayaraghavan  and  Keith  (1990(c)) 
for  the  modified  algorithm.  The  approximate  factorization 
technique  was  found  to  be  robust  and  efficient.  The  unique 
advantage  of  these  techniques  is  that  w.'Ji  the  same  unsleaJy 
formulation,  both  time  accurate  unsteady  results  and  fast 
convergent  asymptotic  steady  state  solutions  can  be 
obtained  in  less  computer  time  than  by  using  a  steady  state 
formulation 

Conclusion 

Recognition  of  the  mathematical  similarities  between 
internal  transonic  flow  and  cavitatcd  bearing  flow  is  an 
important  and  useful  concept.  To  the  best  of  authors’ 
knowledge,  such  similarities  have  not  been  pointed  out 
before.  The  analogy  has  been  exploited  to  determine  the 
characteristics  within  the  cavitatcd  region  and  the  jump 
conditions  across  any  discontinuity  in  the  flow  field.  By 
virtue  of  similarities  between  the  two  flows,  advanced 
concepts  of  transonic  flow  computations  can  be  incorpo¬ 
rated  into  numerical  predictions  of  cavitation  in  bearings. 
Determination  of  the  reformats  rn  boundary  using  both 
shock  fitting  and  shock  capturing  methods  arc  discussed. 
With  the  conservative  formulation  of  the  governing  ‘uni¬ 
versal  equation’,  the  shock  capturing  method  is  very 
effective  and  simple  to  implement.  The  concepts  of  tran¬ 
sonic  flow  computation  have  been  developed  and 
successfully  incorporated  in  three  important  areas,  namely, 
algorithm  development,  grid  arrangement  and  control, 
and  efficient  solution  computation.  The  results  obtained 
arc  encouraging  and  it  is  believed  that  many  more  such 
extensions  may  be  possible  which  will  result  in  improved 
numerical  prediction  of  cavitatio  i  in  bearings. 
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